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Ice sheets appeared in the northern hemisphere around 3 million years ago and 
glacial-interglacial cycles have paced Earth's climate since then. Superimposed on 
these long glacial cycles comes an intricate pattern of millennial and sub-millennial 
variability, including Dansgaard-Oeschger and Heinrich events. 

There are numerous theories about theses oscillations. Here, we review a number 
of them in order to draw a parallel between climatic concepts and dynamical system 
concepts, including, in particular, the relaxation oscillator, excitability, slow-fast 
dynamics and homoclinic orbits. 

Namely, almost all theories of ice ages reviewed here feature a phenomenon of 
synchronisation between internal climate dynamics and the astronomical forcing. 
However, these theories differ in their bifurcation structure and this has an effect 
on the way the ice age phenomenon could grow 3 million years ago. All theories 
on rapid events reviewed here rely on the concept of a limit cycle, which may be 
excited by changes in the ocean surface freshwater balance. The article also reviews 
basic effects of stochastic fluctuations on these models, including the phenomenon 
of phase dispersion, shortening of the limit cycle and stochastic resonance. It con- 
cludes with a more personal statement about the potential for inference with simple 
stochastic dynamical systems in palaeoclimate science. 

Keywords: palaeoclimates, dynamical systems, limit cycle, ice ages, 
Dansgaard-Oeschger events 

1. Introduction 

The Pliocene and the Pleistocene cover approximately the past five million years. 
The climatic fluctuations that characterized this period may be reconstructed from 
numerous natural archives, including marine, continental and ice core records. These 
archives show a complex climate history. Ice sheets appeared in the northern hemi- 
sphere around 3 to 3.5 million years ago [TJ [2]. The volume of these ice sheets 
fluctuated with the variations of the seasonal and spatial distributions of incom- 
ing solar radiation (insolation), which are induced by changes in the geometry of 
the Earth's orbit and the angle (obliquity) between Earth's equator and the eclip- 
tic [3HS] . This is called the astronomical forcing [f] Glacial cycles had an average 

f The astronomical forcing will generally be taken into account here in the form of a normalised 
measure of insolation during the month or on the day of summer solstice at a northerly latitude, 
typically 60 or 65° N. This is a fairly complex, aperiodic signal, with dominant harmonics corre- 
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Figure 1. Climatic fluctuations over the late Pleistocene. Ice ages are reconstructed using 
the oxygen isotopic ratio of the calcite shells of benthic foraminifera [9] . Within the last ice 
ages, large temperature fluctuations were recorded in Greenland, here depicted by changes 
in the oxygen isotopic ratio of ice [10| 

duration of about 40,000 years [BJ until about 800,000 years ago. The dominant 
period of glacial cycles increased around 800,000 years ago and this is referred to 
as the Middle Pleistocene Transition. Data and models about the Middle Pleis- 
tocene Transition are reviewed in ref. [7J. Time-series analyses based on band-pass 
filtering provide further evidence of the non-linear nature of the climate response 
to the astronomical forcing, from about 1.4 Myr ago [5). The latest four glacial 
cycles, in particular, are distinguished by a pronounced saw-tooth time-structure: 
ice accumulates over the continents during about 80,000 years and then melts in 
about 10,000 years (Figure [TJ. 

Superimposed on these long glacial cycles comes a complex pattern of millennial 
and sub-millennial variability |llj . For example, the Greenland record features at 
least 20 events of abrupt rise and slower decline in oxygen isotopic ratio (a proxy 
for temperature) O [13] and methane [14] during the latest glacial epoch. These 
events are known as Dansgaard-Oeschger events. They were found to occur from 
at least the last glacial inception [TS] and the Antarctic ice core record provide 
evidence that they are characteristic of Pleistocene glacial climates [IB]. Some of 

sponding to the phenomena of precession (23716, 22428 and 18976 years); and obliquity (41000 
years) [5] 
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these events follow pulses of iceberg discharges into the North Atlantic Ocean, 
called 'Heinrich events' p~7l - H9j . Heinrich events and Dansgaard-Oeschger events 
have left climatic footprints all over the globe |20| . including in Antarctica [16] . The 
current interglacial period is referred to as the Holocene. It is also characterised by 
millennial and centennial variability, mainly observed in the North Atantic |21H23j , 
but of a much weaker amplitude than during the preceding glacial period. 

The present paper reviews attempts to explain these fluctuations with concepts 
that originate in dynamical system theory. These are the concepts of limit cycle, 
synchronisation and excitability. The central message of the paper is that current 
theories of ice ages and rapid events may often be interpreted in terms of generic 
deterministic models, which are also used in other areas of Science like biology 
and ecology. However, stochastic parameterisations are an essential part of any 
complex system model, and their effects on climatic oscillations have to be taken 
into account. 

Dynamical system theory entered palaeoclimate science with idealised models 
representing the response of ice sheets to the astronomical forcing. These models 
were directly derived from the physics of the ice-sheet-atmosphere system [241/12 7j . 
Ghil and Childress [21], in particular, insisted on the interest of analysing such 
models in terms of bifurcation theory. For modelling the complex carbon cycle re- 
sponse authors sometimes adopted a more heuristic approach by considering simple 
models and confronting the results to palaeoclimate evidence [2H] . 

Nowadays climate research is largely oriented towards large climate simulators 
(typically: general circulation models) , which are developed to include as many cli- 
mate processes as possible. However, thinking in terms of dynamical system theory 
remains insightful. Indeed, the behaviour of a complex system at a certain spatio- 
temporal scale is in practice often dominated by a few leading modes, of which the 
dynamics may be captured fairly convincingly with a low-order dynamical system. 
Climate scientists are increasingly using this property. For example, they formulate 
simpler models to explain the seemingly complex behaviours observed in ocean- 
atmosphere simulators. Examples have been provided in the recent years focusing 
on interannual }30| . centennial |31) and millennial 32, 33 variability. In parallel, so- 
called hysteresis experiments, which aim at identifying the number of stable states 
in individual components of the climate system such as the ocean circulation [34] 
or ice sheets [35] contribute to a dynamical-system founded understanding of the 
climate system. This approach may also help us to predict and communicate about 
the proximity of bifurcations, which may result in catastrophic climatic changes. 
Timmermann and Jin [36] termed predictability of the third kind our ability to an- 
ticipate bifurcation phenomena, by reference to the predictabilities of the first and 
second kind originally introduced by Lorenz [37] , 

The article is structured as follows. Section 2 reviews some of the basic con- 
cepts of oscillator theory. This is no substitute for proper textbook reading, but 
the reader will find essential notions and definitions needed to understand the re- 
mainder. Section 3 reviews how these concepts enter theories of ice ages and rapid 
events. Section 4 discusses effects of stochastic fluctuations and, finally, section 5 is 
a more personal statement about the potential for inference with simple stochastic 
dynamical systems in palaeoclimate science. 
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Figure 2. Sketch of several forms of relaxation oscillations. (A.) The vector space is struc- 
tured by a slow manifold with several stable branches. All points of the state space are 
attracted towards the stable branches of the slow manifold (in full line) along the fast 
direction, which is here the horizontal. The slow evolution consists in an upward or down- 
ward course along the slow manifold depending on whether the system lives on the right- 
or left-hand side of the null-cline of the slow- variable. The relaxation oscillation consists 
in alternate jumps between the two branches of the slow manifold. (B.) Trajectories are 
rapidly attracted towards a region of the phase-space influenced by a saddle point. In the 
scenario dispayed here, there exists a combination of parameters for which the limit cycle 
crosses the saddle point. In that case, the period of the orbit, which includes the saddle-n- 
ode, is infinitely long. It is a 'homoclinic orbit', hence this particular bifurcation is named 
a homoclinic bifurcation. The homoclinic orbit only exists at the bifurcation point, but 
it influences the orbit when the parameter is close to the bifurcation. This is the reason 
why one refers to the 'ghost' of the homoclinic orbit. There is another scenario for which 
different saddle-nodes are connected to each other. The orbit and the associated bifurca- 
tion are then said to be 'heteroclinic'. (C.) The relaxation oscillation is organised around 
a fixed point, with complex eigenvalues with a positive real part. The bifurcation giving 
rise to this orbit is a Hopf bifurcation. (D.) One example of excitability, here depicted for 
a slow-fast system. The system resides in a stable space, but a fluctuation may cause an 
ejection out of the unstable (dashed) branch of the slow manifold. The system then loops 
all the way through the slow manifold before coming back at rest. 
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2. Vocabulary and elementary notions. 

The reader will find an accessible introduction to dynamical system theory and con- 
cepts in ref. |35j. More formal background on oscillator theory, albeit a bit dated, 
is available in ref. [35]. Bifurcation and oscillator theory is explicitly connected to 
climate theory in ref. (see, in particular, chapter 12) and ref. [30], chapter 
7. Background on synchronisation and an introduction to the phenomenon of ex- 
citability is available in ref. [3T] . Finally, the Scholarpaedia peer-reviewed web-site is 
an increasingly rich and authoritative source of information on dynamical systems. 
Only the notions essential for the present article are summarised here. 

Oscillator: The oscillator is a dynamical system that has a globally attracting 
limit cycle. In more simple terms, it oscillates even in absence of an external 
drive. Here, we are interested in oscillators to describe climate phenomena, 
which involve dissipation of energy. The minimal model for a dissipative os- 
cillator includes two ordinary differential equations, of which at least one is 
non-linear. 

Relaxation oscillator: The relaxation oscillator is a particular kind of oscilla- 
tor featuring an interplay between relaxation dynamics (generally fast) and 
a destabilisation process (generally slow). The relaxation is the process by 
which the system is attracted to a region of the phase space. This evokes the 
relaxation of a spring. In a relaxation oscillator the system continues to evolve 
slowly after the relaxation phase. During this slow evolution phase the system 
stability diminishes gradually until the system is ejected out of its relaxation 
state, either towards another relaxation state, or to the same relaxation state 
via a dissipative loop. In this review we will encounter three kinds of relaxation 
oscillators (Figure [2]) : relaxation founded on slow- fast dynamics (involing a 
slow manifold); relaxations structured by a homoclinic orbit (involving only 
one relaxation state), and relaxations structured around a focus. More details 
are given in the caption of Figure [2j 

Excitability: An excitable system has a globally attracting fixed point (it does 
not oscillate spontaneously). However, an external perturbation may have 
the effect of exciting it. During this excitation, the system is being ejected far 
from its fixed point and then returns to it. 

Link between relaxation dynamics and excitability: In practice it is often 
found that a relaxation oscillator may be transformed into an excitable system 
by a mere change in parameter, and vice-versa. The reason is the following. 
A relaxation oscillation is often structured globally in the phase space, for 
example by a slow manifold (Figure [2]A_) or by one or several saddle points 
Figure [2^3). Suppose now that the oscillation displayed by such a system 
ceases because a parameter has been changed. The system is then no longer 
an oscillator, but the 'backbone' of the oscillation dynamics are still latent 
in the phase space because the elements that structured the limit cycle (the 
slow manifold or the saddle points) have not disappeared. Consequently, the 
system may be run on a trajectory close to the defunct limit cycle if it is being 
pushed by some external force (the excitation) into the region of the phase 
space previously occupied by this limit cycle. This point is illustrated on the 
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basis of slow-fast dynamics on Figure [2j3, but similar excitation dynamics 
generally occur near any kind of 'explosive bifurcation', that is, bifurcations 
that give rise rapidly to a fully developed limit cycle. This includes homo- 
clinic, heteroclinic, and certain Hopf bifurcations (two examples follow and 
are illustrated on Figure [6]). 



3. Oscillators, relaxation and excitability in palaeo climates 

(a) Models of ice ages 

(i) The Saltzman et al. models. 

Saltzman established a theory in which ice ages are interpreted as a limit cy- 
cle synchronised on the astronomical forcing. Saltzman and his collaborators wrote 
a series of articles on the subject, starting with the introduction of the limit cy- 
cle idea [32] and synchronisation hypothesis [43] . the interpretation of the Middle 
Pleistocene Transition as a bifurcation [33], and the more complete models in the 
mid-1990s [e.g. ref. [45]. The full theory is developed in a book [40 . Here, we con- 
centrate on two intermediate models [46j[47]. They are called SM90 and SM91, by 
reference to the authors (Saltzman and Maasch) and the year of publication. The 
variables I, fi and 9 are the continental ice mass, CO2 concentration and deep-ocean 
temperature, respectively. The reader is referred to the original publications for the 
meaning and value of the different parameters. They are not crucial here; it suffices 
to know that they are all positive. 
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In both models the first equation describes the ice mass response to changes 
in CO2 (/z) and the astronomical forcing (Fj(t)) Saltzman adopts the so-called 
Milankovitch view [f] that an increase in insolation causes a decrease in ice mass. 
Increases in CO2 or in ocean temperature have the same effects. 

The other two equations describe the dynamics of CO2 and the response of deep- 
ocean temperature to changes in ice volume. It is further assumed that the mean 
state of climate varied slowly throughout the Pliocene-Pleistocene, in particular 
in response to a 'tectonically-driven' decline in the average concentration in CO2, 

f In fact, this view is introduced by Murphy 48 but it is developed mathematically in the 
'Canon of Insolation' [4] authored by Milankovitch 
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Figure 3. Bifurcation diagrams of the Saltzman and Maasch models (SM90) and (SM91) 
as a function of _F M , here treated as a constant control parameter. Note the difference in 
scales on both axes. Black lines are fixed points. Continental ice mass (7) is shown as a 
function of tectonic forcing F M . The red lines indicate limit cycles, shown as the minimum 
and maximum values of 7 along the limit cycle. Unstable fixed points or limit cycles are 
denoted by dashed lines. Calculations and figures were made using the pseudo arc-length 
continuation software package AUTO |50| . 




time (fraction of 76 ka) time (fraction of 110 ka) 

Figure 4. Periodic and fixed point solutions for SM90 (with — 1.2ka/ppm) and SM9f 
(with Ffi = 0.5ka/ppm), near the sub-critical Hopf bifurcation points. The dynamics of 
SM90 slow down near the unstable fixed point, while the limit cycle of SM9f is much more 
decoupled from the position of the fixed point. 

consistently with an earlier proposal |49j . This tectonically-driven decline is here 
modelled as a slow decrease in the forcing term throughout the Pleistocene. 

Consider the bifurcation diagrams of the SM90 and SM91 models with respect 
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Figure 5. Two histories of ice volume generated with the same models (SM90) and (SM91), 
using the astronomical forcing and a decline scenario for the tectonic forcing of CO2. The 
scenarios were chosen here to evidence the rise and decline of 100,000-year ice ages and are 
not the same as those used by Saltzman. Observe the explosive character of the appearance 
and decline of 100,000-year cycles in SM91, with early and late excitations of the limit 
cycle by the astronomical forcing, and the smoother character of the evolution on SM90. 
Time is running from right to left. 

to F^, assuming no astronomical forcing (Ft = 0) (Figure pj). The systems are then 
said to be free or autonomous. Depending on F^, both models show regimes with 
a stable fixed point, and regimes for which the fixed point is unstable so that the 
system orbits along a limit cycle. 

These considerations led Saltzman to interpret the Middle Pleistocene Transi- 
tion as a bifurcation between a 'quasi-linear' response regime to the astronomical 
forcing (in the fixed-point regime) to a regime of non-linear synchronisation (reso- 
nance) on the astronomical forcing. He concluded that ice ages would occur today 
even in absence of astronomical forcing. The main effect of the astronomical forcing 
is to control the timing of glaciations. 

Saltzman's theory is seductive because it explains in a consistent framework 
several aspects of the Pleistocene climate history, including the change from linear 
to non- linear regime [5], the presence of 100,000 year periodicity in climate records 
[5"T] . the lack of a 400,000-year spectral peak in the ice- volume record (such a peak 
appears in the simple piece- wise linear model devised by Imbrie and Imbrie |52| . 
due to rectification of the precession signal), the synchronisation of deglaciations 
on the astronomical forcing [53| 54], and the occurrence of large climatic transitions 
even when eccentricity which modulates the effect of precession on insolation, is at 
its lowest. 

The difficulty for accepting Saltzman's models as a definitive theory lies in the 
physical interpretation of the CO2 equation. This equation encapsulates all the 
interesting dynamics of the system and it is thus crucial to the theory. Some semi- 
empirical justification for the CO2 equation is given in ref. [H] but the form of 
this equation has undergone some somewhat ad hoc adjustments in SM90. The 
form present in SM91 is again different, with important effects on the bifurcation 
structure, while the authors did not justify this latter change based on physical or 
biegeochemical considerations. 

To better appreciate the stuctural differences between the two models let us 
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return to the bifurcation diagrams. Consider SM90. As the forcing is decreased the 
fixed point gives rise to a locally unstable limit cycle. The system must therefore 
find a stable limit cycle further away from the fixed point, but in this case not much 
further. This stable limit cycle is under the influence of the unstable fixed point, 
and in particular the system slows down when it passes near it (Figure [4]) . This is 
scenario 'C depicted on Figure [2] The limit cycle evolves as the tectonic forcing is 
further decreased, until it shrinks smoothly around a perpetually glaciated state. 

In SM91, the bifurcation induced by the decrease in tectonic forcing is much 
more explosive. The system lands on a stable limit cycle that turns out to be little 
affected by the position of the unstable point. The cycle dynamics do not show clear 
phases of acceleration and the system cannot be regarded as a relaxation system. 
The limit cycle disappears abruptly as the tectonic forcing is further decreased, 
through a phenomenon called a saddle-bifurcation of cycles. The consequences of 
the difference between the bifurcation structures of SM90 and SM91 may be further 
appreciated in the transient experiments shown on Figure [5] 

(ii) Paillard's (1998) ice age model (P98). 

Paillard has been advocating the concept of relaxation for understanding palaeo- 
climate dynamics, both ice ages and the more abrupt events, since the publication 
of a seminal paper [55] in 1994. We return to this article later on, and concentrate 
on another article published in 1998 |56j . in which Paillard introduces a conceptual 
model of ice ages. Ice volume dynamics respond to an ordinary differential equation: 



In this equation, the ice volume x is linearly relaxed to xr with characteristic 
relaxation time tr. This relaxation process is further perturbed by the astronomical 
forcing F(t) with a characteristic time r/. Such a system is said to be hybrid [57] 
because the relaxation equation involves a discrete state variable, here denoted 
y. Its state may be 'deep glacial' (G), 'mild glacial' (g) or 'interglaciaP (£). The 
numerical values of Xr and tr depend on this climate state. Climate states y follow 
a sequence i —> g —> G according to a set of conditions formulated on the level of 
glaciation x and insolation. Namely, the transition g — > G is triggered when the 
forcing F(t) exceeds a certain threshold. Occurrence of G drives climate quickly 
into an interglacial state i because xr(G) and tr{G) are specified in the model to 
be low. 

Paillard is not very specific about the physical meaning of the discrete variable, 
but it accommodates the paradigm that the Atlantic ocean circulation has gone 
through three different states during the latest glacial period: intermediate cir- 
culation, shut-down of the circulation, and modern, deep-sinking circulation. The 
system (P98) features the concept of slow-fast relaxation dynamics. However, this 
is not an oscillator because the shift from g to G is determined by the course of the 
external forcing. The Middle Pleistocene Transition is induced in (P98) in a fashion 
similar to Saltzman, and on the basis of similar physical assumptions (tectonically- 
driven decline in CO2). The drift in climatic conditions induced by tectonics is 
accounted for by a term added to the astronomical forcing. In a later review, Pail- 



dj£ _ x R (y) - 
dt r R (y) 



x F(t) 



(P98) 
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lard |58j further emphasises empirical evidence for the relevance of the relaxation 
concept in the phenomenon of deglaciation. 



(iii) The Gildor Tziperman model. 

Gildor and Tziperman [55] take a moderate step towards higher model com- 
plexity by considering a slightly more explicit representation of atmosphere, ocean, 
sea-ice and land-ice dynamics. Namely, the ocean is divided into 8 boxes, and the 
atmosphere into 4. Sea-ice fraction responds to standard energy balance equations. 
More crucially, land-ice growth is influenced by a somewhat controversial feedback 
between sea-ice and precipitation. The feedback is controversial because it is as- 
sumed that cold climate results in a reduction in ice volume: sea-ice growth causes 
a reduction in precipitation in ice-covered areas and, by this mechanism, almost 
suppresses accumulation of snow on ice sheets. The latter then no longer compen- 
sates for ice ablation and ice volume shrinks. 

A free oscillation arises from the fact that the ice volume thresholds for switch- 
ing sea-ice cover 'on' and 'off' differ. In other words, sea-ice displays a hysteresis 
response to variations in ice volume. This is exactly the principle of the slow-fast 
relaxation oscillator depicted on Figure [2]^ : The curve of equilibrium of sea-ice 
with respect to ice volume is the slow manifold, and ice volume integrates the state 
of sea-ice in time. In turn, this oscillation can be synchronised on the astronomical 
forcing. 

The Gildor-Tziperman model is coupled to a biogeochemical cycle in a compan- 
ion paper [60], but the essential dynamics of the glacial oscillation are unchanged. 
Tziperman et al. |61| further comment on the model and its property of synchronisa- 
tion on the astronomical forcing, and find that its behaviour is essentially reducible 
to a hybrid dynamical system. 



(iv) The Paillard-Parrenin model. 

Paillard and Parrenin propose yet another relaxation model in 2004 (PP04). 
The prognostic variables are ice volume /, the area of the Antarctic continental ice 
sheet A and the atmospheric concentration in CO2 (jj) (a . . . j are parameters): 
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As in the other ice age models, ice volume is a slow variable driven by the astro- 
nomical forcing. It is here coupled to a variable with a similar time scale (ta ~ 7j) 
and a faster one (t m = tj/3). The term H(—D), where H is the Heaviside function, 
represents the ventilation of the Southern ocean. CO2 is released into the atmo- 
sphere when the Southern ocean is ventilated (D < 0), which drives deglaciation. 
Ice then grows slowly, until a Southern ocean ventilation flush sends the system 
back to interglacial conditions. Ocean ventilation is thus the fast process in this 
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Figure 6. Bifurcation diagrams and phase-space trajectories of the free Paillard-Parrenin 
model (PP04) and the van der Pol (VDP) model. Both systems display explosive bifurca- 
tion scenarios, but the details are different. PP04 exhibits a sub-critical Hopf bifurcation, 
while the limit cycle of VDP is explicitly framed by a slow manifold. Phase space trajec- 
tories are drawn near the bifurcation points, that is : g = 0.4 in PP04 and /3 — 0.9 in 
VDP. The Heaviside function in PP04 is approximated as H(x) = atan(500a;)/7r + 0.5 for 
analysis with AUTO. 



model and it is the only non-linear process accounted for. Though, contrary to the 
Gildor-Tziperman model, it does not present a hysteresis behaviour. Consequently, 
the glacial cycles featured by this model cannot be interpreted in terms of shifts 
between the branches of a slow manifold. 

To better understand the dynamics of glacial cycles in this model we consider 
the bifurcation diagram along typical solutions in the phase space for the free (i.e. 
unforced) system (Figure[6]). The parameter g is taken in this example as the control 
parameter, in order to preserve Saltzman's idea that ice age cycles appear as the 
consequence of a slow perturbation of the carbon cycle. As in SM90, PP04 exhibits 
a limit cycle arising from a sub-critical Hopf bifurcation. The dynamics along the 
limit cycle close to the bifurcation point are strongly influenced by the presence of 
the unstable focus. This is the configuration 'C shown in Figure[2p. Depending on 
g, the focus is either on the low-ice- volume side of the limit cycle (i.e.: the system 
spends most of its time with high CO2) or on the high- volume side of the limit 
cycle (i.e.: the system spends most of its time in low CO2). Parrenin and Paillard 
estimate that we are currently in the second configuration. 
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(v) A minimal model of ice ages. 

It has been claimed [6"ilR)5] that any model that has some form of 100,000 year 
internal periodicity could be used to reproduce the course of ice volume over the 
last 800,000 years. Taking the argument at face value, Crucifix [51] used one of 
the simplest possible slow-fast oscillators: the van der Pol oscillator, with minimal 
modifications to account for the astronomical forcing and the asymmetry between 
the phase of ice build-up and melt during the late Pleistocene (a, f3, 7, and r are 
parameters; F(t) is the astronomical forcing): 



The system dynamics are determined by the structure of the slow manifold 
x = y 3 /3 — y. The parameter (3 controls the position of the fixed-point on the 
slow manifold and, consequently, the ratio of times spent by the system in the two 
branches ('glacial' and 'interglacial') of the slow manifold. The ice age curve can 
be captured with some tuning (Figure [7]) , although it is fair to add that a small 
change in parameters may shift the timing of one or several ice age cycles. This 
minimal model was used to challenge intuitive arguments about the predictability 
of ice ages [64] . 



(i) Dansgaard-Oeschger events as relaxation oscillations. 

Welander [55] introduced the concept of relaxation oscillations in the context of 
ocean dynamics. He described a heat-salt oscillator involving exchanges of heat and 
salt within a single oceanic column, coupled to a phenomenon of surface temper- 
ature relaxation. The destabilisation process needed for the relaxation oscillation 
to appear is here related to diffusion between the deep ocean and the mixed layer. 
The system dynamics are further controlled by the mean freshwater flux at the top 
of the ocean column. It determines the transitions from a regime characterised by 
perpetual convection in the oceanic column, to a regime with intermittent convec- 
tion (oscillation), and finally to a regime with no convection [66) . The bifurcations 
between the different regimes bear the character of global bifurcations, with the os- 
cillation period approaching infinity near the bifurcation points (in particular, the 
second bifurcation bears the character of a homoclinic bifurcation). The heat-salt 
oscillator belongs thus to the class 'B' on Figure [2] 

Welander [57J and Winton and Sarachik [55] later introduce the concept of an- 
other kind of relaxation oscillator in the ocean. It involves the meridional structure 
of the ocean thermohaline circulation, and the key non-linear process is the merid- 
ional advection of heat and salt. The oscillations featured by this model are termed 
'deep-decoupling' oscillations [55J. Given that the slow process now relates to heat 
accumulating in the global ocean, the characteristic time of deep-decoupling oscil- 
lations is of the order of 1,000 years. The net flux of freshwater delivered to the 
North Atlantic acts as a bifurcation parameter controlling the transition between 
non-oscillating and oscillating regimes in the Winton-Sarachik model |69j . 




(VDP) 



(b) Models for millennial climate variability 
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Figure 7. Astronomical forcing, x and y trajectories obtained using system ( VDP I with 
a = 30, /3 — 0.75, 7 = 0.4 and r = 36 ka (1 ka = 1,000 years). Blue dots are an authoritative 
natural archive thought to mainly represent fluctuations in ice volume and deep ocean 
temperature, and compiled in ref. [9]. 



Millennial oscillations have since been observed across a hierarchy of ocean mod- 
els, including 3-D ocean models with prescribed freshwater flux and restoring con- 
ditions to surface temperature |70) , and 3-D models coupled to a simple atmosphere 
[7T1 172] . Sakai and Peltier [73] proposed that millennial deep- decoupling oscillations 
could explain Dansgaard-Oeschger events. Colin de Verdiere et al. [33J[7U[75] com- 
plement this early proposal with a fairly complete theory based on ocean circulation 
model experiments. The oscillations described by Colin de Verdiere et al. involve 
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the processes of turbulent vertical mixing (neglected in Winton and Sarachik |68j ) . 
advection, and convection, which unify the salt oscillator with the deep-decoupling 
oscillation model. Incidentally, Colin de Verdiere [74] dismisses the non-linearity of 
the equation of state as the cause of the oscillations. 

There is, across the model hierarchy, consistency about the fact that the tran- 
sition between the oscillating circulation regime and the so-called diffusive, ha- 
line regime (without deep convection) is associated with a homoclinic bifurcation 
[33l[76]. The nature of the bifurcation between the convective regime and the oscilla- 
tion is more model-dependent. Timmermann et al. [76], based on experiments with 
the 8-ocean-box Gildor-Tziperman model, find a Hopf bifurcation; salt-conserving 
experiments with a 2-D ocean model show a transition towards a finite-period cycle, 
but of increasing period as the bifurcation is approached; experiments with a more 
idealised model, formulated as a 2-equation dynamical system, reveal the signature 
of an infinite-period bifurcation [33 [f] The latter implies that Dansgaard-Oeschger 
events, at the time when they appear soon after the glacial inception process, should 
be very long but of a similar amplitude as the Dansgaard-Oeschger events coming 
later in the glacial cycle. This feature is consistent with the Greenland ice core 
record (Figure [I]). More specifically, the first Dansgaard-Oeschger cycles that ap- 
peared at the beginning of the glacial era were characterised by a long 'plateau' 
phase (also called: interstadial) during which the thermohaline circulation was cer- 
tainly very active [15] ■ In the Colin de Verdiere et al. theory, the plateau phase is 
the phase of the trajectory influenced by the 'ghost of the saddle point' [73]. 

(ii) Dansgaard-Oeschger cycles as the manifestation of an excitable system. 

Given the explosive nature of the bifurcations involved in ocean dynamics it 
is no surprise to find excitability properties in ocean models. Weaver and Hughes 
|70j discuss this effect in salt-conserving experiments with an idealised-geometry, 
ocean model. The ocean-atmosphere model of intermediate complexity CLIMBER 
(CLIMate BiosphERe mode) was shown to exhibit excitability properties when 
boundary conditions are set to be typical of the latest glacial era [77] . The ocean 
circulation has then one stable state, with moderate Atlantic overturning, and a 
'quasi-stable state' with more intense overturning. The conceptual sketch of the 
excitation cycles shown by Ganopolski and Rahmstorf |78] on their Figure 1 can 
be interpreted in terms of slow-fast dynamics, in which the different states of the 
ocean circulation constitute the different branches of a slow manifold. The intense 
overturning state, which is the 'plateau' phase of the Dansgaard-Oeschger event, 
may thus be viewed as the repelling branch of the slow manifold in the excitable 
regime (Figure [2p). The excitable Dansgaard-Oeschger hypothesis was used as 
a possible basis to explain how a weak forcing, exogeneous to the system, could 
explain the observed 1500-yr periodicity of Dansgaard-Oeschger cycles (on this 
periodicity: see ref. [751110] but see the other view in ref. [HI])- Two such theories 
were developed on the basis of experiments with CLIMBER. One suggests that 

f This particular case was not illustrated on Figure [5] There is no saddle point along or near 
the orbit, but there is a combination of parameters for which a fixed point appears on the limit 
cycle. Beyond this particular parameter value, this fixed point splits into a saddle and a node. This 
particular parameter value correpsonds to an 'infinite-period' bifurcation. In practice, as long as 
the limit cycle exists, the trajectory slows down near the point where the saddlc-nodc will appear. 
Some authors then refer to the influence of the 'ghost' of the saddle point [38] - 
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Dansgaard-Oeschger events are excited by stochastic fluctuations, modulated by 
a weak, hypothetical solar periodic forcing [82] (more on the effects of stochastic 
fluctuations in section [I]). The alternative theory suggests that the excitation is 
induced by the interference between two solar forcings with periods close to 1470/7 
(= 210) and 1470/17 (« 87) years [53], possibly combined with noise [54] . 

(hi) Heinrich cycles as a relaxation oscillation. 

MacAyeal [85 proposed an ice-binge/purge theory to explain Heinrich events. 
The theory rests on experiments with a 1-spatial direction model of ice flow dy- 
namics. Suppose, as a starting point, that ice volume grows in response to net accu- 
mulation of snow. The growth continues until the accumulated effect of geothermal 
heat flux causes basal sliding. A volume of ice is then released into the ocean (this 
is the 'purge'), causing the release of icebergs characteristic of Heinrich events. Ice 
volume thus decreases, until ice accumulation wins over so that ice volume can grow 
again. The ice-binge/purge model is thus a relaxation oscillator combining a slow 
integrating process (ice mass accumulation) with a fast lateral discharge process. 

(iv) Coupling between Heinrich and Dansgaard-Oeschger events. 

To what extent Heinrich events may interfere with Dansgaard-Oeschger dy- 
namics? Paillard [SSJ [SB] investigated this question by coupling the MacAyeal ice 
model — but reduced to ordinary differential equations by Galerkin truncation — 
with a 3-box ocean model. The coupling simply assumes that ice released into the 
ocean causes a net freshening of the surface of the North Atlantic that alters the 
deep-ocean circulation. Paillard realised that this coupling could lead to fairly non- 
intuitive and complex effects, such as the succession of Dansgaard-Oeschger events 
of decreasing amplitude between Heinrich events. This succession is known in the 
litterature on palaeoclimate records as Bond cycles [18) . Paillard also found that 
the oscillations are aperiodic in this model under certain parameter configurations. 

The issue is further explored in [69] , based on the Winton-Sarachik ocean model, 
and in [76] , based on the slightly more sophisticated Gildor-Tziperman ocean model 
[55] . The objective was to study the response of deep-decoupling ocean oscillations 
to prescribed Heinrich cycles. Schulz et al. [55] noted that deep-decoupling os- 
cillations could be synchronised on the Heinrich cycles. Timmermann et al. [76] 
then proposed, on the basis of numerical experiments with a fairly idealised model, 
that Heinrich events excite Dansgaard-Oeschger cycles because the variation in ice 
volume caused by a Heinrich event modifies slowly the amount of net freshwater 
released in the ocean. In turn, they suggested, Dansgaard-Oeschger may have a 
control on ice volume growth. This yields a two-way coupling between Dansgaard- 
Oeschger and Heinrich events. 

Experiments with more comprehensive models of the ocean-ice-sheet-atmosphere 
system [571 H5J generally support the idea that the different water and heat fluxes 
involved in the different phases of ice build-up and iceberg release are quantitatively 
sufficient to support a coupling between ice sheets and ocean circulation during the 
latest glacial era. However, it was also noted that " three-dimensional thermome- 
chanical ice-sheet models are unable to satisfactorily reproduce the binge-purge 
mechanism without an ad hoc basal parameterisation." [89j . 
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To address this difficulty a theory in which Dansgaard-Oeschger events trigger 
Heinrich events was recently proposed [89]. The ice shelve plays a key role, in 
blocking the ice stream flow from the ice sheet to the oceans. Heinrich events occur 
when this ice shelve is broken, for example under the influence of ocean sub-surface 
warming associated with a Dansgaard-Oeschger event. The resulting model is a 
system displaying a slow ice-build-up - Heinrich release cycle excited by fluctuations 
in ocean sub-surface temperature. 

(v) Holocene oscillations and relationship with Dansgaard-Oeschger events. 

The much smaller ocean oscillations that characterised the Holocene period 
may also be a relaxation phenomenon. Schulz et al. [31] observe oscillations in the 
atmosphere-ocean model of intermediate complexity ECBILT-CLIO. These oscilla- 
tions are related to the convective activity in the Labrador Seas. 

Schulz et al. [5D] considered the existence of such an oscillator in an earlier refer- 
ence and speculated on the possible interactions between the centennial oscillations, 
millennial oscillations, and Heinrich cycles. They considered a model in which each 
of these three kinds of oscillations is modelled as a Morris-Lecar relaxation oscilla- 
tor. Their working hypothesis is that glacial conditions induce a coupling between 
these oscillators. They then observed that a very stable 1500-yr oscillation appears, 
which they interpreted as a model equivalent of Dansgaard-Oeschger events. 

4. Stochastic effects 

The myriad of chaotic motions that characterise the dynamics of the ocean and the 
atmosphere may be taken into account in the form of parameterisations involving 
stochastic time-processes. The method was introduced in climatology in the 1970s 
[91] and the theoretical justifications, which allow one to model chaotic motions as 
a (linear) stochastic process, are reviewed in ref. [92], [93] . In a statistical inferential 
framework, the stochastic parameterisations may also be viewed as a way to account 
for the distance necessarily existing between the concepts and dynamics featured 
by the model, and the complex system being observed. 

The effects of stochastic processes on relaxation oscillators and excitable systems 
are generally well documented in the literature because this is a topic of general 
interest [M] . Here we review some of them in the specific context of palaeoclimate 
dynamics. 

(a) Stochastic effects on ice age dynamics 

(i) Phase dispersion. 

One of the basic effects of noise on oscillators is the phenomenon of phase disper- 
sion: A weak stochastic forcing on an oscillator causes a fading out of the memory 
of the exact initial conditions, even though the gross structure of the oscillation 
visualised in the phase space is conserved. The phenomenon is well known and it is 
an immediate consequence of the neutral stability of the phase of a free oscillator 
with respect to fluctuations. It was early suggested that this phenomenon of phase 
dispersion may concern ice ages [95) , but it is more commonly believed that ice ages 
are phase-locked on the astronomical forcing. This phase-locking should act against 
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dispersion and permit a very long predictability horizon of ice ages. Though, a phe- 
nomenon of phase dispersion may happen in oscillators that are locked on a periodic 
forcing. A stochastic fluctuation may momentarily cause a burst of desynchronisa- 
tion, called phase slip, during which the system is unhooked from its corresponding 
deterministic trajectory and attracted to another trajectory, which leads or lags 
the original one by one forcing period (ref. [STj . sect. 3.1.3). The difference between 
phase diffusion in a free system and in a periodic-forcing-driven oscillator is that 
the diffusion effect has, in the latter, a quantum nature. More formally, it is said 
that the stochastic forcing disperses the system states around the different attrac- 
tors that are compatible with the forcing. In a work in preparation we suggest that 
the astronomically-forced climate system may satisfy the conditions for a similar 
phenomenon of phase dispersion to occur (B. De Saedeleer, M. Crucifix and S. 
Wieczorek, unpublished data, 2011). Given that the astronomical forcing is aperi- 
odic the description of the phenomenon requires a suitable theoretical framework, 
which relies on the notion of a 'local pullback attractor'. The equivalent of a phase 
slip is, in the aperiodic forcing context, a stochastic shift from one of the determin- 
istic pullback attractors to another one. The phenomenon is illustrated based on 
experiments with the VDP model on Figure [8] 

(ii) Reduction of period. 

Additive fluctuations generally reduce the period of relaxation oscillators. In 
an oscillator presenting a homoclinic orbit such as the Duffing oscillator, additive 
fluctuations reduce the time spent near the unstable focus This implies that 
even at the corresponding bifurcation point in the deterministic system, the return 
time of oscillations in the stochastically-perturbed system remains finite. In a slow- 
fast oscillator such as the van der Pol oscillator, additive fluctuations generally 
result in early escapes of the branch of the slow-manifold on which the system lies 
(Figure[9]). The period of the oscillator is thus affected by a correction that increases 
approximately linearly with the noise variance in the slow-fast van der Pol oscillator 
97 [f] This property was used at least once in Pleistocene theory, in the silicate 
weathering hypothesis advanced by Toggweiler [05]. Additive fluctuations reduce 
the limit-cycle period from 800,000 years to about 100,000 years. The reasons for 
the period reduction being so dramatic are left for another article. 

(b) Stochastic effects on Dansgaard-Oeschger dynamics 

(i) Stochastic excitation and resonance. 

Noise may naturally act as an excitation agent in an excitable system (Figure [9]). 
The topic is extensively reviewed in ref. |94j . Excitation loops are sporadic if the 
noise amplitude is weak, in which case the recurrence time of the events is set by 
the noise amplitude, while the amplitude of the events is set by the structure of 
the deterministic vector field. The frequency of excitation loops increases with the 
noise amplitude, until the system behaviour is qualitatively similar to a limit cycle 

f This result is established assuming fluctuations added to the slow variable. A much more 
general theory, suitable for different slow manifolds and additive fluctuations to slow and fast 
variables is now available |98| . 
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Figure 8. The phenomenon of trajectory shift illustrated with the model (VDP) with 
astronomical forcing, with parameters r = 36 ka, j3 — 0.75, 7 = 0.4, a = 30. Black: the 
deterministic system. The generated history reproduces reasonably well the fluctuations 
of ice volume of the last 800,000 years. Red: same system but with an additive Wiener 
process added to the fast variable y, with variance b = 0.2(ka) _1//2 . One possible realisation 
of the stochastic system is shown. Observe the solution shift around 450,000 years ago. 
Depending on the realisation the shift may occur at different places. 

regime. This is the coherent resonance regime. For yet higher noise amplitudes, the 
limit cycle structure is destroyed. 

The concept of stochastic excitation has been considered several times in Dans- 
gaard-Oeschger theories. The idea is introduced based on experiments with an ocean 
general circulation model with idealised geometry and forcing [70]. The effect of 
stochastic fluctuations is not only to excite oscillations in a system normally at rest, 
but also to reduce the period of these oscillations when the system is in oscillatory 
regime. 
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A. noise on excitable slow-fast system 




Figure 9. Two possible effects of additive noisy fluctuations in a slow-fast system. (A.) 
If the slow-fast system is in the excitable regime, noisy fluctuations can cause excitations. 
Excitations will be sporadic if the amplitude of the oscillations is weak (as shown here), 
or regular, as in a limit cycle, in the coherent resonance regime. (B.) If the slow-fast 
system is oscillating, additive fluctuations may cause early or delayed escapes from the 
slow branches, compared to the deterministic system. The net effet is a shortening of the 
cycle duration. 

A phenomenon of non-autonomous stochastic resonance may occur if the noise is 
superimposed to a weak external drive. For this to happen the autonomous system 
needs to be stable but excitable. The external forcing must be too weak to cause 
excitation by itself. The role of the noise is to provide the additional power to induce 
excitation. The timing of the excitation is then related to the phase of the external 
forcing. The mechanism was proposed several times [7H [55] to explain the 1500- yr 
recurrence time of Dansgaard-Oeschger events. The idea remains questioned, either 
on the ground that the 1500-yr recurrence time observed in palaeoclimate records is 
coincidental |81j or on the ground that the 1500-yr external forcing is unidentified 
[33} [83] . A more subtle case of stochastic resosonance involves the combination 
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of noise with two solar cycles of 210 and 87 years, which yields the concept of 
'ghost resonance' |100j for which some support, albeit not conclusive, is found in 
the observations [80] . 

(ii) Decreased sensitivity to noise in resonant oscillators. 

Coupled oscillators may exhibit, collectively, a resonance period that is more 
robust to external fluctuations than the uncoupled oscillators. Schulz et al. [5U] 
used this property to explain the stability of the Dansgaard-Oeschger recurrence 
period of 1500 years in presence of random fluctuations, without having to invoke 
an external forcing. 

(hi) Pseudo-oscillations in two-well systems. 

Finally, a behaviour reminiscent of oscillations may occur in a system that is 
neither oscillating nor excitable, but which presents several stable states. Noise then 
simply induces jumps between these different states. The simplest mathematical 
model is the Langevin equation and this is on this basis that Schulz et al. [31] 
interpret the Holocene oscillations observed in the ECBILT-CLIO climate model. 

5. Concluding discussion : Can dynamical systems be used 

for inference? 

The review has shown that relaxation oscillations are a popular and powerful model 
to explain oscillations observed in the Pleistocene record. The concept of relaxation 
implies some form of slow-fast separation, in the sense that at least one component 
of the system spends most of its time in 'quasi-equilibrium' states (this may be a 
'slow manifold branch' or a region influenced by a saddle-point, depending on the 
system structure), with acceleration phases. 

Some of these models were constructed following a fairly careful procedure of 
truncation of a system of partial differential equations, which describes some of 
the fluid dynamics of the climate system. Others were proposed on a more con- 
ceptual basis, the idea being precisely to test a hypothesis based on palaeoclimate 
observations. The latter approach is sometimes criticised, on the ground that box- 
models, for example, cannot reasonably be taken as an adequate representation of 
the complex dynamics of the oceans [101 j . 

This leads us to the last question of this review: can dynamical systems be 
used for inference on palaeoclimates? Inference implies that something is being 
learned by confronting a model to observations. This inference process may take 
the form of a calibration procedure (update our knowledge on parameters on the 
basis of observations) or a model selection procedure (which model, among different 
alternatives, explains the observations best). 

The position taken here is that there is not such a thing as an 'attractor' of 
the climate system that is to be 'discovered'. The hope is that some of its modes 
of behaviour are sufficiently de-coupled from the rest of the variability to justify 
the fact simple dynamical systems may capture the fundamental dynamical prop- 
erties of these modes, and we want to learn about these modes from palaeoclimate 
observations. 
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The programme is challenging. Indeed, it was underlined that different physical 
assumptions may lead to dynamical systems with dynamical properties that are 
similar enough to produce a convincing visual fit on palaeoclimate data |61j . The 
message is largely echoed in the present review. The modeller's challenge is there- 
fore to operate a model selection on more stringent criteria than just fitting some 
standard time series. For example, palaeoclimate observations may yield constraints 
on the bifurcation structure of the system. The Middle-Pleistocene Transition is an 
attractive test case in this respect. 

In a statistical inference process, the observations should be a plausible outcome 
or realisation of the model. This makes sense only if the model has a stochastic com- 
ponent, which describes its uncertainties, limitations, and the noise that emerges 
from the chaotic motions of the atmosphere and oceans. 

Stochastic dynamical systems begin to be used for inference on palaeoclimate 
time series. In a method called 'potential analysis', the climate system is modelled 
as a Langevin equation, that is, the combination of a down-gradient drift with 
a Wiener additive process, and inference is made on the number of wells of the 
potential function [102J. The method was applied on Pleistocene climate records, 
yielding the conclusion that the number of wells increased from 2 to 3 over the 
course of the Pleistocene |103j . 

However, our position so far has been to favour a Bayesian methodology, because 
it allows one to encode physical constraints in the form of prior distributions on 
model parameters. The Bayesian formalism is also naturally designed for model 
calibration, selection, and probabilistic predictions. 

The fact is that Bayesian methods for selection and calibration of dynamical 
systems on noisy observations are only emerging. In a recent attempt we considered 
a particle filter for parameter and state estimation [104] . To be honest, there is 
ample room for progress. Whether the process of inference with simple dynamical 
systems on palaeoclimate data will lead new insight in this context still needs to 
be demonstrated. 
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